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Abstract 


WAHLS, RICHARD A. Development of a Defect Stream Function, Law of the 
Wall/Wake Method for Compressible Turbulent Boundary Layers (Under the direction of 
Dr. Fred R. DeJamette) 

The method presented is designed to improve the accuracy and computational 
efficiency of existing numerical methods for the solution of flows with compressible 
turbulent boundary layers. A compressible defect stream function formulation of the 
governing equations assuming an arbitrary turbulence model is derived. This formulation 
is advantageous because it has a constrained zero-order approximation with respect to the 
wall shear stress and the tangential momentum equation has a first integral. Previous 
problems with this type of formulation near the wall are eliminated by using empirically 
based analytic expressions to define the flow near the wall. The van Driest law of the wall 
for velocity and the modified Crocco temperature-velocity relationship are used. The 
associated compressible law of the wake is determined and it extends the valid range of the 
analytic expressions beyond the logarithmic region of the boundary layer. The need for an 
inner-region eddy viscosity model is completely avoided. The near-wall analytic 
expressions are patched to numerically computed outer region solutions at a point 
determined during the computation. A new boundary condition on the normal derivative of 
the tangential velocity at the surface is presented; this condition replaces the no-slip 
condition and enables numerical integration to the surface with a relatively coarse grid using 
only an outer region turbulence model. The method has been evaluated for incompressible 
and compressible equilibrium flows and has been implemented into an existing Navier- 
Stokes code using the assumption of local equilibrium flow with respect to the patching. 
The method has proven to be accurate and efficient. 
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1 INTRODUCTION 


It is generally accepted that turbulence is the most complicated kind of fluid motion. 

Many researchers have directed extensive effort towards the physical understanding and 

computation of turbulent flows. Turbulence has been described by Hinze 1 as "...an 

irregular condition of flow in which the various quantities show a random variation with 

time and space coordinates so that statistically distinct average values can be discerned." 

The major problem when attempting to compute a turbulent flowfield from first principles 

is that the time and space scales of the turbulent motion are extremely small. The 

computational grid required to fully resolve such a flow in this manner are beyond the 

2 

limits of todays computer technology. Anderson et al provide an estimate of the spacing 
required for a typical flowfield in which "10 5 points may be required to resolve just 1 cm^ 
of the flowfield." It is because "statistically distinct average values can be discerned" that 
modeling techniques can be developed that allow us to solve turbulent flowfields for 
engineering purposes. 

Currently there are two basic approaches for computing turbulent flowfields . The 
first approach is known as a large eddy simulation. This approach attempts to compute the 
maximum amount of information about the turbulent motion as it attempts to capture the 
turbulence on a real-time, real- space basis. Large eddies, which are responsible for the 
majority of the momentum transport, are resolved on the scale of the grid. Subgrid 
modeling is used to describe the effects of the small eddies, which cannot be resolved with 
practical grids. It has been projected that fully converged large eddy simulations for 
simple airfoils is beyond the capability of current computers. 

The second, and more prevalent, approach is to calculate the mean motion of the 
fluid by use of the time-averaged Navier-Stokes equations. These equations are often 
referred to as the Reynolds equations. They are derived by decomposing the dependent 


variables into mean and fluctuating components and then averaging the equations over time. 

The resulting equations are expressed primarily in terms of mean variables with the 

transient effects of turbulence being described by the fluctuating quantities that appear in the 

turbulent or Reynolds stress terms; the Reynolds stress terms require modeling. This 

approach does not attempt to capture the turbulence on a real-time, real-space basis but 

instead relies on a turbulence model to convey the effects of the turbulent motion on the 

mean flow. Also, the solution schemes do not need to be time accurate and are able to 

converge much more rapidly than with large eddy simulations. On the other hand, 

information on the dynamics of turbulence is lost. It is fortunate that the simpler statistical 

4 

approach is adequate for general engineering purposes . 

Turbulence modeling is the most important and most difficult aspect of the statistical 
methods. The turbulence model is necessary to mathematically close the system of 
governing equations and must fully describe all of the effects of the turbulent motion on the 
mean flow. Several excellent review papers on the state of turbulence modeling have been 
given by Rubesin , Lakshminarayana , Marvin , and Rubesin and Viegas . Turbulence 

3 

models can be categorized as either zonal or global in reference to their range validity . The 
zonal model is one which is designed for a specific application and is usually developed 
from this specific class of experiments. In general, these models do an excellent job of 
predicting flows of the particular class but begin to lose accuracy and break down when 
extending beyond these flows. The global models, however, involve more complex 
functions, such as the field equations for the turbulence quantities, and have a wider range 
of applicability. The more complex models require more empirically based coefficients 
than do the zonal models; these coefficients are determined from more than one class of 
experiments. The basic necessity of empirical coefficients restricts the global models from 
successfully predicting flows of all types. At present, a true global, or universal, model 
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which is accurate for all flows has yet to be determined. Thus, it can be stated that current 
turbulence models are restricted to various zones of application. 

Turbulence models can also be classified as either eddy-viscosity or stress- transport 

6 8 
models . Eddy-viscosity models rely on the Boussinesq concept, which models the 

turbulent shear stress as the product of an effective viscosity and a mean rate of strain. The 

effective viscosity is proportional to the velocity and length scales of the turbulence in the 

various regions of the flow. Eddy-viscosity models can be further classified according to 

the method of determining the effective viscosity; in other words, the various eddy- 

viscosity models are functions of the method used to describe the velocity and length scales 

of turbulence. They are classified as either zero-, one-, or two-equation models, where the 

designation refers to number of field equations used in addition to the usual mean flow 

equations. The simplest is, of course, the zero-equation model, which determines the eddy 

viscosity based entirely on the properties of the mean flow. The two primary examples of 

this type are those of Cebeci and Smith and Baldwin and Lomax . These models are 

widely used in practical engineering applications for simple shear flows and in many 

Navier-Stokes codes 5 . The one-equation models typically rely on the solution of the 

turbulent kinetic energy equation to represent the velocity scale, and the two-equation 

models incorporate the solution of an additional equation which represents the length scale 

of the problem. The one- and two-equation models are more general than the zero-equation 

models but, in many instances, do not improve the accuracy of the solutions. 

The stress-transport models use the Reynolds-stress equations to model turbulence 
in the mean flow equations. Here closure is achieved through the solution of the mean 
turbulent field and requires the solution of up to six additional differential equations for the 
Reynolds stresses, thus significantly increasing the computational effort required beyond 
the capability of current computers 5 . As a result, several attempts to simplify the stress- 
transport equations have been made and are discussed in detail in reference 5. These 
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approaches make assumptions, some arbitrary and others based on the physics of the flow, 
that allow the stress-transport equations to be reduced to an algebraic form. Note that even 
in differential form, these equations still require modeling of several terms which involve a 
number of empirical coefficients. The main advantage relative to eddy-viscosity models is 
that the stresses are able to respond immediately to changes in the rate of strain. These 
models are appropriate for a wider range of flow conditions including separated and 
recirculation regions for which eddy-viscosity models are not valid. 

The current investigation is concerned with the solution of compressible turbulent 
boundary layers over solid surfaces. The flowfields studied are restricted to attached 
boundary layers under adiabatic wall conditions. The formulation is derived such that it is 
consistent with the eddy-viscosity concept. Over the years, considerable research has been 
devoted to the development of numerical methods for just such situations. The 
computational effort required of these methods is highly dependent on the number of grid 
points necessary to accurately resolve the flowfield. Boundary-layer flows in general 
require a finely spaced grid in order to accurately calculate their properties and their 
influence on the external flow. Turbulent boundary layers, in particular, demand a very 
fine grid near the wall in order to resolve the high gradients of the physical properties in 
this region; however, a much coarser grid similar to that necessary for a laminar boundary 
layer can be used in the outer region away from the wall. Blottner 1 1 presented an excellent 
example demonstrating the different grid requirements for laminar and turbulent boundary 
layers. Blottner's figure is reproduced here in figure 1. With a common outer edge of the 
boundary layer defined and the requirement of one percent accuracy of the shear stress at 
the wall, Blottner estimated the number of uniform intervals necessary to resolve laminar 
and turbulent boundary layers on a flat plate. He found that the turbulent profile required 
twenty times the number of intervals as the laminar profile. The additional intervals 
required are a direct result of the high gradient region near the wall, and to a lesser extent. 
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the increased thickness of the turbulent boundary layer. An obvious solution to this 
problem is the use of a nonuniform grid which clusters grid points where the most 
resolution is needed; in this case, grid points are clustered in the high gradient region near 
the wall. Nonuniform grids are in everyday use and have been for many years. Despite 
this technique, the grid requirements still become excessive when dealing with complex 
two-dimensional and even the simplest three-dimensional flows as computational time can 
take many hours or days even on the fastest available computers. 

An additional technique which has become popular recently is the use of wall 
functions. In effect, wall functions replace the highly clustered grid points and 
corresponding numerical computations near the surface with empirically based, analytic 
expressions for velocity, temperature, and other quantities as necessary for compatibility 
with a particular turbulence model. It is, in fact, the potential reduction of computational 
effort that is the main incentive behind the development of wall function techniques. 
Experimental investigations of turbulent flows, and the corresponding analytic 
investigations, have established expressions for the near-wall region which have proven to 
be quite accurate and robust. The use of such expressions should not be a concern as it is 
apparent that all turbulent flow analyses are based on empiricism to some extent. The 
logarithmic law of the wall for velocity is one of the most widely used and accepted 
expressions of this type. It is, in fact, one of the most widely accepted empiricisms in all 
of fluid mechanics and has proven useful even beyond the bounds of its strict assumptions. 

The basic concept of the wall function methods is to determine analytically the 
flowfield variables from the surface to a point near the surface where the analytic 
description is valid. This is the first grid point off of the surface. From this point outward, 
the flowfield is resolved with an appropriate numerical method and associated turbulence 
model. The analytic description of the flowfield variables at the first grid point provides 
boundary conditions used to compute the remainder of the flowfield. Most methods simply 
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patch the analytic and numeric solutions at this point while some, such as the recent method 
of Walker, Ece, and Werle , formally match the outer limit of the inner solution to the 
inner limit of the outer solution. A common characteristic of current wall function 
approaches is the precomputation designation of the first grid point off of the wall; that is, 
the section of the flowfield to be determined analytically is decided prior to the 
computation. This characteristic is independent of the numerical method and turbulence 
model. The grid point location is crucial because it must be within the range of validity of 
the analytic model at all streamwise locations. This grid point is typically chosen to lie 
within the viscous sublayer or logarithmic region of the boundary layer. As such, this 
approach does not take full advantage of the analytic expressions because their use is 
limited to only a portion of their valid range. It follows that near-wall turbulence models 
are still required although not to the same extent as without the use of wall functions. 

3 

Rubesin's review of turbulence modeling points out several major advantages that 
wall function methods have over methods that integrate to the surface. First of all, wall 
function methods have proven to be quite economical due to the reduction in the number of 
grid points required to resolve the flowfield as well as to the increase in the allowable time 
step due to the increased size of the minimum grid spacing. Secondly, present near-wall 
turbulence models have inherent physical uncertainties and, in some instances, result in a 
numerically stiff set of equations near the surface. The use of the wall functions has been 
shown to increase the accuracy of the solutions near the surface in many cases while 
relieving this near-wall stiffness. As far as the accuracy is concerned, it should be 
advantageous to model the inner region with an empirical representation which can be 
validated by direct experimental observation rather than with another model, such as an 
eddy-viscosity model, which cannot be validated by direct experimental observation. The 
point is that error can enter through a model if said model is based on indirect empiricism 
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such as velocity measurements being correlated and put in the form of an eddy viscosity; 
the eddy viscosity is not a property which can be measured directly. 

One of the most effective wall function methods to date is that of Viegas, Rubesin, 
and Horstman , which is an extension of the method originally developed by Viegas and 
Rubesin 14 to increase computational efficiency for their study of shock wave/boundary 
layer interactions. Their method has proven to be successful not only by increasing the 
efficiency of their solution procedure by an order of magnitude with respect to the computer 
time per solution, but also by improving the accuracy of their solutions. The method has 
proven successful over a variety of flowfields including two-dimensional separated flows. 
However, several parameters, such as the skin-friction coefficient and shock locations, 
remain somewhat sensitive to the choice of the first grid point off of the surface (the patch 
point of the analytic and numeric solutions). More recently Wilcox 15 reports the sensitivity 
of the location of the first grid point on the skin friction. 

The method put forth in this investigation combines traditional wall function ideas 
with features designed to further improve computational efficiency and address some of the 
shortcomings of previous methods. The improvement in efficiency will be a direct result of 
further grid reduction, and the primary shortcoming addressed is the choice of the location 
of the first grid point. The method proposed considers the turbulent boundary layer as a 
composite consisting of an inner region near the surface with an outer region beyond. The 
two regions are similar to those encountered with zero-equation turbulence models, which 
use separate eddy-viscosity models in the different regions. The inner region is to be 
resolved in a completely analytic manner while the outer region is resolved numerically. 

The inner region exploits the law of the wall for velocity in conjunction with an associated 
law of the wake. The law of the wake enables the effects of streamwise pressure gradients 
to influence the region near the wall. The energy equation is replaced by a modified form 
of the Crocco temperature-velocity relationship. The outer region is formulated and 
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resolved in terms of the defect stream function. This formulation properly characterizes the 
outer portion of the boundary layer as it relates the local mean velocity to the velocity at the 
outer edge of the boundary layer. 

As stated earlier, the law of the wall for velocity is probably the most widely 
accepted empiricism in fluid mechanics. The most common form of this law is known as 
the logarithmic law of the wall and is valid in the fully turbulent, logarithmic region of the 
inner part of the boundary layer. Others have devised laws that are valid from the wall 
through the viscous sublayer and the logarithmic region. One such example for 
incompressible flow is the law given by Liakopolous 1 **, which was used in a prior 
investigation 17 ' 18 by the present author. Far and away the most popular compressible law 
of the wall is that described as the effective, or generalized, velocity approach of van 
Driest 19 , which effectively relates the velocity in a compressible fluid to a corresponding 
velocity in an incompressible fluid. The van Driest law is used in the present method and is 
discussed further in section 2.2. 

Law of the wake empiricisms have not been developed or used to the extent as law 

of the wall empiricisms, particularly for compressible flows. The most popular law of the 

wake is that presented by Coles 20 for incompressible flow. Moses 21 also gives a 

commonly used incompressible law of the wake. Several other laws for incompressible 

flow are given by White in reference 22. McQuaid 23 proposed an incompressible law that 

used a modified shear-stress velocity as the normalizing variable, such that pressure 

gradients from the most favorable to separation could be handled with relative ease, while 

the more common laws of the wake simply use the shear-stress velocity as the normalizing 

variable. The most popular form of a compressible law of the wake appears to be the 

incorporation of Coles' law into the van Driest effective velocity approach. Maise and 

McDonald 24 merely speculated on this extension without a rigorous derivation and found 

25 

this to be a valid approach. This extension was used in the work of Alber and Coats and 
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further scrutinized by Mathews, Childs, and Paynter also with favorable results. 

Squire used a compressible flow extension of McQuaid's law and found the law to give 
questionable results with Mach numbers of 3.6 and above. 

It was felt that an additional contribution of the present investigation could be a 
compressible law of the wake determined directly from compressible flow solutions rather 
than relying on correlations with a law that is based entirely on incompressible flow data. 
This has been accomplished and is discussed in detail in the upcoming sections. The 
precise law of the wake determined here is valid through the inner region and into the inner 
portion of the outer region. The significant effect of including a law of the wake is the 
extension of the region of the boundary layer that is described by the analytic expressions 
beyond the logarithmic region. This provides an additional grid reduction leading to 
increased computational efficiency. Physically speaking, the law of the wake allows for 
the influence of the streamwise pressure gradient on the inner region of the boundary layer. 

no OQ 

The importance of this effect has been discussed by many including McDonald , Patel , 
the present author in references 17 and 18, and more recently by Wilcox in reference 15. 

A primary feature of the present method is the approach used to patch the inner and 
outer solutions. Specifically, the location of the patching is not predetermined by the 
researcher, but rather is determined as part of the computation. This approach has several 
advantages. First of all, the arbitrary choice of the location of the the patching is 
eliminated. Secondly, the choice is made as part of the computation in a manner which 
extends the use of the analytic expressions at each streamwise location. Various means of 
implementing such an approach are discussed in section 3.3. 

The formulation in the outer region is also a key factor with the present technique. 

A prior investigation by the present author ’ formulated the outer region in terms of the 
conventional stream function; the outer region velocity was nondimensionalized by the 
velocity at the edge of the boundary layer. The critical implementation feature was the 
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algorithm developed to interact the inner and outer region treatments. Solutions were 
obtained for incompressible flows over flat plates, and two-dimensional ellipses and 
circular cylinders. Once pressure gradient effects were included in the inner region 
treatment, successful computations were made, although a fair amount of numerical 
difficulty was encountered. The present formulation reduces the numerical difficulty by 
replacing, when possible, computational steps in the interaction routine with equivalent 
analytic steps, This process was facilitated by reformulating the outer region in terms of 
the defect stream function; the outer region is characterized by a velocity deficit, relative to 
the edge velocity, normalized by the shear-stress velocity. 

Much of the outer region analysis to follow is based on previous analyses by 
Clauser 30,31 and by Mellor and Gibson 32 . Clauser's analysis resulted from a desire to 
study turbulent boundary layers excluding upstream history effects. The comparison is 
made to the self-similar laminar solutions of Blasius, and Falkner and Skan, as shown in 
most fluids textbooks, where the governing equations reduce exacdy to ordinary 
differential equation form, resulting in a solution independent of the streamwise location. 
The problem for turbulent flows is much more difficult, as not only must the mean velocity 
profiles be similar, but also the profiles of the turbulence quantities. Clauser observed that 
experimental tangential velocity profiles for incompressible turbulent boundary layers on a 
flat plate reduce to a single curve, independent of the Reynolds number, when plotted in 
terms of the velocity defect normalized by the shear-stress velocity. Several of Clauser's 
figures from his 1956 paper 31 are reprinted here in figure 2 for the purpose of 
demonstration. Figure 2a shows several turbulent velocity profiles from flat plate 
experiments plotted in the conventional velocity ratio format. Figure 2b demonstrates the 
collapse of the data to a single curve when plotted in the velocity defect format. Clauser 
further proved that velocity defect profiles of turbulent boundary layers with streamwise 
pressure gradients are also self-similar if the pressure and skin friction forces are, in his 
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terminology, in equilibrium. Clauser determined, with much effort, that a pressure 
gradient parameter normalized by the ratio of the (incompressible) displacement thickness 
to the shear stress at the wall designates the various families of similar solutions for 
turbulent boundary layer flow. A similar parameter for compressible flow is determined in 
the present investigation. 

Another of Clauser's contributions was the definition of a highly accurate eddy- 
viscosity model valid in the outer region. The model is a result of Clauser's observation of 
the laminar-like behavior of the outer region of the turbulent boundary layer. To 
demonstrate this similarity, Clauser solved the Blasius equation with various nonzero wall 
velocities, as shown in figure 2c, and compared to the turbulent profiles shown in figure 
2b. The similarity away from the wall, as shown in figure 2d, is remarkable. Clauser 
determined the proper velocity and length scales to be the edge velocity and the 
(incompressible) displacement thickness, respectively, and a corresponding constant which 
completed an outer region eddy-viscosity model that remains constant normal to the 
surface. These scales were subsequently verified for compressible flow by Maise and 
McDonald 24 and variations of the model are still widely used today. 

Clauser did not derive the governing equations in terms of the defect stream 
function, but rather used the conventional stream function. Although others studied 
equilibrium turbulent boundary layers after Clauser, it was not until ten years after Clauser 
that Mellor and Gibson put forth a clear, accurate derivation of the governing equations in 
terms of the defect stream function variables. It should be noted that a significant portion 
of the work described in reference 32 was originally presented four years earlier in 
Gibson's doctoral dissertation . The work of reference 32 is restricted to the analysis of 
incompressible, equilibrium turbulent boundary layers. Mellor and Gibson obtained an 
extremely accurate, approximate solution to the governing equations in the limit of 
vanishing shear-stress velocity to edge velocity ratio. They also determined a completely 
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analytic solution of the governing equations by using this approximate equation and taking 
advantage of the first-integral property of the tangential momentum equation. With this 
solution, they established the extremes of Clauser’s pressure gradient parameter from the 
most favorable pressure gradient to separation. The two major problems with their 
formulation concerned the enforcement of law of the wall behavior near the wall and the 
implementation of the wall-layer eddy-viscosity model. Also, the wall boundary condition, 
which is inversely proportional to the shear-stress velocity to edge velocity ratio (the ratio is 
a small number), presents numerical difficulties. The computational complexity Mellor and 
Gibson encountered near the wall probably explains the general lack of popularity of the 
defect stream function formulation. Even Mellor abandoned this formulation for his 
analysis of nonequilibrium, incompressible turbulent boundary layers. The present 
technique completely overcomes the difficulty near the wall by use of analytic expressions 
in the inner region. 

A significant accomplishment of the present investigation is the extension of the 
defect stream function formulation to nonequilibrium, compressible turbulent boundary 
layers. The derivation follows in chapter 2. Chapter 3 presents the inner region treatment 
and discusses its implementation. Finally, chapter 4 discusses the results of the present 
technique. The equilibrium class of boundary layers has been used extensively over the 
years in the development and testing of analytic and numerical methods. The process is 
continued in the current investigation: section 4.1 shows solutions for equilibrium, 
incompressible flow; and section 4.2 discusses equilibrium, compressible flow solutions. 
Section 4.3 discusses the initial implementation of the method into an existing Navier- 
Stokes code and presents results for several cases of compressible flow over a flat plate. 
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2 DEFECT STREAM FUNCTION 
FORMULATION 

32 

The basic formulation parallels that of Mellor and Gibson for incompressible, 
equilibrium, turbulent boundary-layer flow. However, the present treatment is for 
compressible, nonequilibrium, turbulent boundary layers; and the eddy viscosity model for 
the present treatment is much more general than that of reference 32. In fact, it is 
unnecessary to specify the eddy viscosity in the inner layer at all. 


2.1 Basic Equations 

The Reynolds-averaged continuity and tangential momentum equations for 
compressible turbulent boundary-layer flow in two dimensions are 

d(pu) | d(pv) _ q 

dx dy 


du du 

pu— +pv— 
dx dy 



dr 

dy 


( 2 . 1 ) 


where x and y are the tangential and normal coordinates, u and v are the respective velocity 
components, p is the density, and the subscript 'e' designates values at the edge of the 
boundary layer. The shear stress T is 

t = p ~~ - ptPV (2.2) 

dy 

where (1 is the molecular viscosity of the fluid and (-puV) is the Reynolds stress term, a 
direct result of fluctuations inherent to turbulent flow. The Reynolds stress term must be 
modeled because an exact, analytic form is not known. In the present investigation the 

O 

Boussinesq eddy viscosity concept is used rather than a more complex stress transport 
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equation model. The Boussinesq concept imitates the Newtonian assumption that the shear 
stress is linearly proportional to the mean rate of strain; the Reynolds stress is modeled as 

-puv=n,^ (2.3) 

dy 

where p t is the eddy viscosity. The eddy viscosity is proportional to a velocity scale and a 
length scale; the choice of these scales distinguishes the various eddy viscosity models 
from one another. The present derivation models the sum of p and p t as 

p + Pt = K(x,y)pu e 5* (2.4) 

where K is a general nondimensional function of x and y and 8* is the incompressible 

displacement thickness. The incompressible displacement thickness 8! is sometimes 

35 

referred to as the velocity thickness and is defined as 


5t=J(l -£>dy 

31 

This model (equation (2.4)) is of the form proposed by Clauser for the outer region of 
incompressible turbulent boundary layers. The generality of the present form is achieved 
by the use of the arbitrary function K(x,y), which allows this model to duplicate other 
existing turbulence models without affecting the present derivation. Although the function 
K(x,y) can have separate definitions in the inner and outer regions of the boundary layer, 
only an outer region definition is required with the present technique. The present solutions 
were computed assuming K(x,y) = k in the interest of simplicity, where k is the Clauser 
constant. 

The energy equation is accounted for by use of the modified Crocco temperature- 

36 

velocity relationship. The basic relationship was derived by Crocco with following 
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assumptions: steady flow of a perfect gas, Prandtl number (Pr) of one, which implies a 
perfect balance of viscous dissipation and heat conduction, streamwise pressure gradient of 
zero, and a constant value of the specific heat coefficient Cp. The basic relationship is 

T = T w + (T aw -T w )i-^ 


where T w is the wall temperature and T aw is the adiabatic wall temperature. Baronti and 

Libby have made a detailed investigation of the accuracy of the Crocco relation for 

adiabatic flows and state that they found deviations of less than plus or minus four percent 

in the static temperature ratios obtained throughout the boundary layer. Gran, Lewis, and 
38 

Kubota have compared this relation to nonadiabatic experimental cases and found no 

22 

significant deviation. White shows a modification to this equation in which a recovery 
factor is introduced as follows: 

T = T w + CTaw-T w )^-^ 


where r is the recovery factor and is defined as r=Pr^ for turbulent flow (r=Pr^ for 
laminar flow). This equation has proven to be a very good approximation' of the energy 
equation even beyond the bounds of the strict assumptions. 

Adiabatic wall conditions are assumed in the present study; the modified Crocco 
temperature-velocity relationship for adiabatic walls is 



w 


l + r^V 


(2.5) 


where y is the ratio of specific heats (taken to be 1.4), a is the speed of sound, and M is the 
Mach number. The Prandtl number is assumed to be 0.72 in this study. The simplicity 
and effectiveness of this relationship avoids the use of the differential energy equation and 
the corresponding need for temperature laws of the wall and wake. 
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30 31 

In this treatment, the defect stream function of Clauser ’ is modified to account 
for compressibility. The defect stream function f(£,T[(x,y)) of this formulation is defined 

such that 


u 


(2.6) 


where the transformed coordinates are 

y 


^ = x 



The coordinate transformation incorporates a density weighting integral similar to that 
originally proposed by Mager . This particular transformation of the normal coordinate 
was chosen because it reduces to the form used by Clauser, and then Mellor and Gibson, 
for incompressible flow. The shear-stress velocity u* is defined as 



where p w and t w are the density and shear stress at the wall. The compressibility 
transformation yields the density-weighted velocity thickness 8*, which appears throughout 
the derivation and is defined as 


oo 



The boundary layer defect thickness parameter A is defined as 
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oo 


oo 



Au* = u e 5* 

Partial derivatives with respect to x and y are of the form 

_c) 5 _ p 3 

3x 3£ 9x 3r| 8y Ap e 3r| 

2.2 Law of the Wall and Wake 

It is assumed that a law of the wall and wake for velocity in the inner part of the 
boundary layer is known. This law is of the form 

-T = g(y + ,Me) + hCp.rj.Me) (2.7) 

u 


where g is the law of the wall and h is the law of the wake. The inner variable y + is 
defined as 


+ = u2 I= Rcs 1 
v„ “ 1 


where 0) and T) are defined as 

8| _ y 

co = — ■ and ri = ■*- 

K A 


and the Reynolds number based on the edge velocity, the incompressible displacement 
thickness, and kinematic viscosity at the wall is defined as 


u e 8j „ „ 
Res* = co 


* * A 

u A 


V w 
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The compressible pressure gradient parameter, which reduces to Clauser's pressure 
gradient parameter for incompressible flow, is defined as 



where p is the pressure. A modified form of the parameter P will be defined later and 
shown to be the compressible equilibrium parameter. For nonequilibrium flow, this new 
parameter is a function of x and hence for equilibrium flow, this parameter is constant. 

In the present treatment, van Driest's^ effective velocity (or generalized velocity) 
approach, as simplified for adiabatic flows and modified by the recovery factor, is used for 
the law of the wall: 

g(y + .Me) = A P=“(^9 sin { a/ ^ (£) ( ^ y + + B) ) (2.8) 

\ rCr-1) u K 

where a aw is the adiabatic wall speed of sound and k and B are the law of the wall 
constants. Van Driest's derivation parallels Prandtl's 40 derivation for incompressible flow. 
The assumptions used are that the shear stress in the fluid is constant and approximately 
equal to the shear stress at the wall and the mixing length is linearly proportional to the 
distance from the wall. Van Driest incorporated a variable density into the derivation in the 
form of the Crocco relationship. Again, equation (2.8) has been simplified for flows with 
adiabatic wall conditions. Because the ratio (u*/a a w) is generally small, the law of the wall 

is often used as 

g(y + .Me)= i lny + + B + 0[(£ ; ) 2 ] 

1C 


Neither of these expressions is valid in the laminar sublayer. 
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Law of the wake empiricisms have not been developed and used to the extent that 
law of the wall empiricisms have, particularly for compressible flow. It is widely accepted 
that the leading term of the incompressible law of the wake is proportional to y 2 . There is, 
however, disagreement over the order of the second term . Moses , for example, uses a 
second term proportional to y^ while Coles’^ second term is proportional to y 4 . If the law 
of the wake needs to be valid from the wall through the inner part of the outer region of the 
boundary layer only, a single term proportional to y 2 is sufficient to define the wake 
function. In the present compressible treatment, a functional form proportional to y 2 for 
the law of the wake is assumed and is modeled after the leading term of the Moses 
incompressible law of the wake. The present law of the wake for compressible flow is 

h(P,fj,M e ) = — w q 2 (2.9) 

K 

where w = w(fJ,Me) is a coefficient determined as part of the computation. The procedure 
for determining this coefficient is discussed later, and an analytic expression is given as a 
result of the present computations. 

2.3 Shear-Stress Velocity Ratio 

The nondimensional shear-stress velocity is defined as 



In the present treatment, the ratio y is evaluated with the laws of the wall and wake. From 
equations (2.6) and (2.7) it is seen that the velocity u at the match point T| m between the 
inner and outer regions of the boundary layer can be expressed as 

u m = u*[ g(yp + h(7W ) = u e + u*f (^m) ( 2 -10) 

The equation for y is obtained from equation (2.10) as 
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Y_ [ gm + h m - fm ] (2<11) 

In the analysis that follows an expression for the gradient of y with respect to B, is needed. 
Differentiation of the above expression yields 

Y=-7 2 [gm + h m -i r A,] (2.12) 


where 


() = 




The derivative of the law of the wall (equation (2.8)) with respect to % is 


e=-[ i +!t+t-— ]+o[(f-n 


* o 

U v2_ 


K y °e A V 


w 


a aw 


(2.13) 


Substitution of equation (2.13) into (2.12) yields 

■ -A 

Ue Y_ K [ 1 i Uc ( ^ Vw y J + Pe A (hm ~ fm) 

uJd+Z) ie A ’ v » P » P (1 + 1) 

K K 


31 

This term was neglected by Clauser but was found to be important and included by 
Mellor and Gibson in their incompressible analysis. The present compressible analysis 
will also require the following relationships: 


Pe 

Pe 



v w Pw Pw 
V w Pw Pw 
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|iw 


Pw 


1+3 fp~ 
[ 

2(1 + =f-) 



Pw , (r - 1) I Ue 
Pe 1 ^ 


£^=-[m2+2(i-£^) 

Pw Pe 


1 + 


Pe. 

Pw 


)] 


Ue 

Ue 


These relationships are derived from the perfect gas equation of state, the modified Crocco 
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relation for adiabatic wall conditions (equation(2.5)), and Sutherland's viscosity law 

Ue Vw 

where s = 199°R. Note that the term is a known function of Me. 


2.4 Governing Equations 

Now one equation for the defect stream function f can be written. The gradients of 
u with respect to x and y are 

^■= ^ (1 + Yf’ ) + ue jfyf’ ) + Uetf ' ~ 
dx UA dx dx 

9u _ pUeTf " 

3y p e A 

respectively. From the definition of the stream function 

d*¥ J 3^ 

pu = — and -pv = — 

dy dx 


it can be shown that 
= p e u e A (q + 70 
and the flux pv is obtained as 
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(2.14) 


. 3 3n 

-pv = ^(peUeA) (T) + yf) + p e UeA — (yf) + p c u e A (1 +yf’):p- 


In the defect stream function formulation, the tangential momentum equation, equation 
(2.1), is 


® P£_ [ K (£-) 2 r]’ + ^2.(2 f + yf 2 ) - (1 + — -) [(n + yf)f" 

P Pw Pe Pw u e A 

<-) (h . 

+ — - (f +Yf ,2 -7ff")] + M|(ri+70f" +— — — — (f +yf ,2 -yff") 


(1 + 1 ) 

K 

<?> 






(1 + V W 

K 


(2.15) 


= ^I[3L t ^5L. r 3f] 

Pwp 8s 3s 3s 


where s is the nondimensional tangential coordinate defined as 


\ 


s = 



2.5 Boundary Conditions 

There are two surface boundary conditions and one far-field boundary condition. 
From equation (2.14) for the flux pv, the normal flow boundary condition at the surface, 
v = 0 at T) = 0, is observed to correspond to 

f(s,0) = 0 
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The far-field boundary condition is obtained from the definition of the boundary-layer 
defect thickness A: 


A = - 


‘p 3f . . f 3 f . 

— — dy = -A I —dll = 


fa* 

Pe dn 


d 


dll = A (f(s,0) - foo(s)) 


It follows that 

foo (s)= -1 (2.16) 

This boundary condition requires that u approach u e as y approaches infinity. The final 
boundary condition involves the shear stress at the wall. The shear-stress boundary 
condition is 

.. . . 3u 

t w = hnuj (|i + pt) — 

This boundary condition can be written as 

p w „»*- 1^. KfcDP 

T|-*0 p e A 
or 

Mt Kft.11) (— ) J r = (2. 17) 

Tl->0 p e 0) p e 

This boundary condition replaces the usual no-slip boundary condition at the surface 
(u(s,0)=0) and states simply that the shear stress must approach the wall shear stress as y 
approaches zero. 

It should be noted that these boundary conditions reduce to those used by Mellor 
32 

and Gibson for incompressible flow. 
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2.6 First Integral of Governing Equation 

It is always desired to reduce the governing equations to their simplest form 
analytically before turning to numerical analysis. The advantage of such an analytic 
integration is twofold: the order of the equation to be solved numerically can be reduced, 
and boundary conditions can be absorbed and not imposed numerically. The defect stream 
function governing equation has a first integral, and one boundary condition is absorbed 
analytically. As will be seen, however, a limited advantage is gained for the full equation. 
It is only upon making the zero-order approximation in the limit of vanishing shear-stress 
velocity ratio that the full advantage is gained. 

Equation (2.15) can be written as 


[-£S-K(P-) 2 f" 

P pw Pe 


<T> 


+2££-f-(i+^) {nr - (1 ' Yf )f } + — — >f 


Pw 


u e 


(1 +-) (1 +-) Pw 
K K 


+ Me V + j-]' 

Y. K ’ V w P w p 0S pw P 

u e 


(i + ^> 

K 


Pw (5 3s Pw (3 3s 




(1 + 1 ) 

K 


Ue 


pw K 


u e 


t Ue ^w Pe \ 0 7 Pe /dhrn df m \ i Pe 7 


.2 


Using the boundary condition given in equation (2. 17), the above equation can be 

Ue A 

integrated across the boundary layer to evaluate the quantity 1 h . This value is 

• A 

u e 
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l+-T=-[( 1 + 1 )(r + — (2-lG)-Mt(l-iO) — ^< 2 -|G- 


K p p w 


(l4l)Ue Vw 


, . 1 


= 1 +■ 


where the defect shape factor G is defined as 


G = J f 2 dr| 


It is necessary to incorporate a Mach number scaling effect into the coordinates in 
order to determine the conditions for compressible equilibrium flow. The following 
transformation allows for this effect: 


s = ( — ds and rj = 

fl Pe 


— r] 

Pe 


Applying this transformation and integrating across the boundary layer, the first integral of 
the governing equation for arbitrary r\ is 


(l-yf ) §■= n(£-f Kf" - p (1 + i- M?) (rjf - 


d A \ 

S pe 


<nJ) 

K 


— P(l-7f)f(M|-^) + 2p££-f-l +— Y — [p(l+l-M| + -^.)(l-l) 


( 1 + 1 ) 


(1+ 1 ) 

K 


Pw K 


- 2 — p (M| - — — -i^-) + 27 (^p-- A /“■“£“)] ff 

K v w p w ds \ p e ds 0 
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A 

■n 


A 


—7(1- + Y f f 2 dfy ] - YTrlf 2 

Pw Pe Pe 0 dso 

, Y ds \ pe ds 


(1 +“) 
K 


where the prime denotes partial differentiation with respect to T]. The term (1 + — ) is 
defined in equation (2, 1 8), and the parameter y is defined as 




Y 


Note that for the adiabatic wall conditions of this investigation, the density ratio (p w /pe) has 
values between zero and one, which assures that Y is less than or equal to Y- Expressions 
for the ratios <n and (pe/p) can also be written in defect stream function form: 


00 



For edge Mach numbers from incompressible to supersonic this parameter is small so that 
the ratios (§T/Sv) and (Pe/p) are approximately one. However, because e has the limiting 
form 


26 



limit e = 

>oo 1 L 


yMe 

, 4 1 

it is not small for large values of the edge Mach number. The term e is similar to Coles' 
parameter CfM^> which remains finite for large Mach numbers. The skin friction coefficient 
Cf is equal to (2t w /pooUoo 2 ). 

It is apparent that the analytic integration shown in equation (2.20) provides an 
advantage over equation (2.15), the equation before integration. The governing equation, 
which is now in the form of an integro-differential equation, has been reduced to second 
order (leading derivative of f") and one boundary condition (equation (2.17)) has been 
absorbed analytically. The integral terms in equation (2.20) are of higher order, and are 
evaluated easily when solving this equation. The full advantage of the first integral 
equation is seen in the next section where the zero-order approximation is made. 

2.7 Zero-Order Approximation 

As mentioned previously, the shear-stress velocity ratio y is generally small. As a 

result, the dependent variable of the defect stream function formulation f can be expanded 
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in terms of the shear-stress velocity ratio. Mellor and Gibson used this expansion in their 
incompressible analyses to obtain a zero-order, asymptotic form of the governing equations 
with respect to the shear-stress velocity ratio. The expansion of f is 

f = f o + yfi +y 2 f 2 + •••• 

Asymptotic forms are obtained for equations (2.18) and (2.20). For hypersonic flow, 
these asymptotic forms are not the strict zero-order forms because the term e appears in 
both co and (pe/p) and is not negligible. The constrained zero-order forms of equations 
(2.18) and (2.20) are 
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( 2 . 21 ) 


T3T- ®(-) 2 Kfo + (1+2P) T|fj ■ f 0 - 1 
as p e 

A 

where the pressure gradient parameter {3 is defined as 

a-^) 

p,£«.p[i._^«-(i-£w ( i.r))] 

Pw Zr Pe 

Note that equation (2.21) is linear from incompressible to supersonic flow because 0) and 
(pe/p) are approximately one. An additional transformation is necessary to obtain a 
governing equation which is linear in the hypersonic range. This transformation is 
described in reference 42 and shown in Appendix A. 

The full advantage of analytically integrating the governing equation is now 
apparent. The terms in equation (2.20) are either simplified or eliminated when making the 
zero-order approximation; the integral terms of equation (2.20) are considered to be small 
and are neglected. The result is simply a second order differential equation and one 
boundary condition has been absorbed analytically. 

2.8 Equilibrium Flow Approximation 
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The equilibrium condition, as defined by Mellor and Gibson , occurs 
mathematically when the profile of f, and hence u, depend only on the nondimensional 
normal coordinate and not on the streamwise coordinate. These authors noted that the 
streamwise partial derivatives are zero when the coefficients of the governing equation are 
independent of the streamwise coordinate, and they showed that this condition can be met 
exacdy for the zero-order incompressible equation and approximately when higher order 
terms are included. This same situation pertains for compressible flow. If K does not 


A n 

depend on s, P is constant, and both to and (pe/p) are approximately one, the coefficients 
of equation (2.21) are independent of s and, since the boundary conditions are also 
independent of s, the derivative on the left side of equation (2.21) is zero. The more 
complicated coefficients of equation (2.20) depend weakly on s. 

The equilibrium form of the full equation (equation(2.20)) is 


(Y) 

co(— ) 2 Kf" - p (1 + M?) (qf* - — P(l-7f’ ) f (m 2 - ^L) 


Pe 


(i4) (i4) 
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u e 
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Pw (1+ T^ m p w K K y e V w p w 0 
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+^Jf2 d ^] =0 
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( 2 . 22 ) 


The equilibrium form of the constrained zero-order equation (equation (2.21)) is 

co(£-) 2 + (1+2(3) Tjf^-f 0 -l=0 (2.23) 

Pe 

It will be shown that numerical solutions to the nonlinear, equilibrium equation (equation 
(2.22)) and solutions to the approximately linear form of equation (2.23) are virtually 
identical. Each of these equilibrium equations is an ordinary differential equation. 

Mellor and Gibson have shown that the zero-order incompressible equation is in 
the form of a confluent hypergeometric equation, which has an analytic solution . They 
showed equilibrium solutions to exist from P = -0.5, the most favorable pressure gradient 
case, to P = oo, the incipient separation case. The corresponding compressible analysis is 
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presented in reference 42, and the analytic solution for compressible flow is given in 
Appendix A. 
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3 Inner Region Treatment 


It is this treatment which most distinguishes the present method of solution from 
previous methods. In effect, the treatment replaces the inner region numerical 
computations with an empirical representation. The inner region generally encompasses the 
inner-most twenty percent of the boundary layer and can extend beyond the logarithmic 
region of the boundary layer. It can also be thought of as the region where an inner eddy- 
viscosity model is used in the zero-equation eddy-viscosity approach. The need for an 
eddy-viscosity model in this region is completely eliminated. The empirical representation 
is in the form of the law of the wall and the law of the wake for velocity as defined in 
equation (2.7) and the modified Crocco relation for the temperature as defined in equation 
(2.5). The laws of the wall and wake need only be valid in the inner region for the present 
method; the Crocco relation is used across the entire boundary layer. The outer region of 
the boundary layer is where the outer eddy-viscosity model pertains. 

In general, there is one point, the "match point," where both the laws of the wall 
and wake and the outer eddy-viscosity model are correct. The present treatment assures 
that the derivatives of the defect stream function through f " are continuous at this point. 

The term "match point" should not be construed to mean that the inner and outer solutions 
are being matched in the formal sense. Actually, these solutions are being patched at one 
point. As discussed in references 17 and 18, this procedure is completely analogous to 
patching the inner and outer eddy viscosity models in the zero-equation modeling approach. 

It should be noted that the inclusion of the law of the wake in the inner layer model 
means that the match point is not confined to the logarithmic part of the boundary layer. 

Several key points must be addressed when deriving and implementing this inner 
region treatment. First of all, an equation relating the velocity, as represented by f ' , and its 
integral, represented by f, through the inner region must be available to insure a continuous 
solution across the match point. Secondly, an equation which allows the determination of 
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the match point as part of the computation must be determined. Finally, a procedure must 
be developed which allows the implementation of the patching of the inner empirical and 
outer numerical solutions in an efficient manner. These points will now be addressed. 

3.1 Equation relating f and F 

An equation can be established which relates f and f' throughout the inner region of 
the boundary layer. This equation follows from the definitions of f and the laws of the wall 
and wake: 


f = 


•i 

/ 


df J df 

— dq =n — - 

dq dq 


q 


(Og_ + dh_) 

dq 


'dg_ 

dq 


dq 


A relationship between the q and q is needed to form the final inner region relationship 
between f and f' . Such a relationship is 


H 

T\ 


— {l+ - 2X< tane - 2X)] + 0[(£) 4 ] } 


Pe 


cos 2 0' 


(3.1) 


where the terms 0 and % result from the modified van Driest law of the wall (equation 
(2.8)) and are defined as 


1 = 



1 . Pw 

Pe 


and 


0 = XK(-lny + + B) 
K 


Note that the term %k is no larger than the shear-stress velocity ratio y and the term 0 is 
typically small, which enables use of the small angle approximation in equation (3.1). The 
relationship between f and f ' is written as 
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The lowest order approximate relationship between f and f ' is written as 


"H 


e (— — i\ C j— .^.ivu* \2i 

f= ’ l( ^T ) - v [( ^ )] 
J * 


(3.2) 


This equation pertains throughout the inner region, where the empirical law of the wall and 
wake is valid. Note that because h is proportional to rp, the integral term is proportional to 
Tp; the integral term is small throughout the inner region because fj m is small. The use of 
this equation at the match point insures the continuity of f and f 1 . The match point is 
positioned so that the derivative f" is continuous. The transformation given in equation 

(2.19) is easily applied to the equations of this section. 

3.2 Match Point Location 

The match point location is determined using the equilibrium form of equation 

(2.20) , which is the first integral of the tangential momentum equation as given in equation 
(2.22), and equations (2.7) and (2.8) for the law of the wall and wake, which are used to 
evaluate the f " term in equation (2.20). The equilibrium flow assumption can be used 
because the derivatives with respect to s in equation (2.20) are small in the inner region 
where the match point is located. Equation (3.2) is used to relate f and f . Equation (2.22) 
is evaluated at the match point and, with some manipulation, the constrained zero-order 
form of the governing equation is 
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(3.3) 


K0) PSL{— i- + (^) m } + (l + 2K<\ m © m -f m -l=0 
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Using the transformation given in equation (2.19) and equation (3.1), the lowest order 
relation between r\ and Tj is 

■n=n(— ) 3/2 

Pe 

Because the law of the wake in equation (2.9) is of order Tj 2 near the wall, the 
integral term in equation (3.2) is of order T|3 and can be neglected as nm is small. If 
equation (3.2) is substituted into equation (3.3), the equation for the match point location is 

ATtf-iW + £<»£a=0 (3.4) 

K Pw 

where the constrained zero-order form of the term A is 

A = (Pw)3/2(I /^ +2 ^(|.) m ) + 12w K(i) Pm (3.5) 

Pe K y pw ^ Pw 


It is in equation (3.4) that the law of the wake plays its most significant role in the analysis 
and, from equation (3.5), it is the leading coefficient of the law of the wake which is truly 
significant. The solution for f) m is easily determined from equation (3.4) as 
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(3.6) 


3.3 Implementation of Inner Region Treatment 

There are several methods with which to implement this inner region treatment. 
The essential feature is that the inner region formulation properly interact with the 
numerically computed outer region solution. As discussed earlier, the formulation is such 
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that f, f' , and f " are continuous at the interface between the two regions. In primitive 
variable terminology, this means that the integral of the velocity profile, the velocity u, and 
the velocity gradient du/dy are continuous across the match point. In effect, the inner 
region formulation provides boundary conditions for the outer numerical solution. This is 
the same concept that has been used in other wall function methods. Again, the main 
differences of the present method are the increased region of the boundary layer that is 
modeled empirically and the fact that the size of this region is not predetermined by the 
researcher (ie. previous methods limit the use of empirical expressions to a predetermined 
distance from the wall) but is determined as part of the solution process. 

Once the match point is determined, empirical expressions are evaluated at this point 
so that the properties of the flow are determined and used as boundary conditions for the 
numerical computation of the outer region. These new boundary conditions replace the 
traditional no-slip at the wall boundary conditions normally enforced at the wall for viscous 
flows. The region of the boundary layer requiring numerical computation on a highly 
clustered grid is significantly reduced. This is true for previous wall function methods and 
was, in fact, the main incentive behind the development of these methods. As mentioned 
before, the present treatment increases the savings. 

At a given streamwise location, the present treatment yields a single, physically 
accurate match point location; this location, however, is variable in the streamwise 
direction. Problems with the implementation of the new boundary conditions arise because 
the grid is typically established prior to the numerical computations; it would be fortuitous 
if the variable location of the match point would always, or ever for that matter, fall on an 
existing grid point. Computationally, there are several options as to the point of application 
of the empirically determined, match point boundary conditions. 

The obvious first option is to apply an adaptive grid technique that forces the grid 
point nearest the wall to coincide with the match point at each streamwise location; the first 
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grid point off of the wall adjusts itself with each iteration, or time step, during the 
computation such that it is always located at the current match point. The match point 
boundary conditions are enforced at the first grid point. This approach is very complex for 
the general case of nonequilibrium flows, but is relatively simple for equilibrium flows. 
This is the primary approach used in the present investigation when computing equilibrium 
flows (both incompressible and compressible); the details of this technique as applied here 
for equilibrium flows are given later. 

An approach must be developed that relates the match point to a general, fixed grid. 
Several options present themselves. Since the inner region of the boundary layer is thin 
relative to the boundary layer thickness, it should be admissible as a first approximation to 
impose the match point boundary conditions at the geometric surface. This is comparable 
to translating the point of application of boundary conditions in thin airfoil theory to the 
airfoil chordline. A more accurate approximation is to apply the match point boundary 
conditions at the grid point closest to the match point. Flowfield information at the grid 
points between the match point and the surface is supplied by the law of the wall and wake. 

Another possibility is the definition of a new transformed normal coordinate that is 
defined such that the match point is located at a constant value of the new coordinate for all 
streamwise locations. Such an approach requires the definition of a scaling function used 
in the transformation of the normal coordinate. The scaling function is dependent on the 
streamwise location and varies from flowfield to flowfield. Melnik 44 and Walker, Ece, and 
Werle have developed methods incorporating this approach. 

The notion of applying surface slip velocity boundary conditions is reminiscent of 
31 

Clauser 1 s analysis of the outer region of turbulent boundary layers and comparison to 
laminar boundary layers with slip velocities. Note the slip velocity corresponding to the 
outer region shown in figure 2d. It would be desirable to use information determined from 
the inner region treatment and the matching process to define surface boundary conditions 
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which enable the proper numerical computation of the outer region of the boundary layer. 

In effect, the outer region is assumed to extend all the way down to the wall so that only an 
outer eddy-viscosity model is required. The shear stress at the wall associated with the 
outer eddy-viscosity model is forced to be equal to the physically accurate wall shear stress. 
As a result, the numerical computations from the wall to the match point yield nonphysical 
information. However, the numerical solution becomes physically accurate upon reaching 
the match point, thus providing the proper description of the outer region. The condition 
specified in this approach is on the velocity gradient at the wall, which implies a condition 
on the shear stress at the wall. In words, the numerically computed wall shear stress must 
be equal to the physically accurate wall shear stress which is determined from the empirical 
inner region and matching treatments. Mathematically, the boundary condition on the 
velocity gradient enforced at the wall is 


/3u\ _ _ p w u 
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(3.7) 


and is written in defect stream function variables as 


r — {—) 1 _ P w /Pe\ 2 St 
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This equation is a direct result of the definition of the shear stress (equation(2.2)) using the 
eddy-viscosity concept (equation(2.3)). The arbitrary function K is chosen in the form of 
an outer region model, such as K=k where k is the Clauser constant. This outer solution 
includes a nonzero surface slip velocity and the associated surface density p s . 

Conditions on the velocity gradient at the wall have been used prior to this 
investigation by other researchers. The condition used by Gorski et al is common among 
wall function approaches. The condition is of the form 
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and is simply the derivative of the common logarithmic law of the wall. This condition 
does not enforce the outer region model at the wall, but rather provides the compatibility 
between the velocity, as determined from the law of the wall, and its normal derivative in 
the near wall region. Melnik 44 , on the other hand, has developed a condition for 
incompressible flow that is also in the tradition of Clauser. His condition, which was 
developed independently from the present investigation, is essentially the same condition as 
the incompressible form of the present boundary condition. 

It is felt that the implementation of the inner region treatment will be best handled in 
the general, nonequilibrium case by use of the slip-velocity approach just described. This 
approach is used in the primitive variable applications of this investigation. It has also been 
successfully implemented into the equilibrium boundary layer solution procedure using 
defect stream function variables. 


4 Results and Discussion 


This chapter traces the development of the present method in chronological order. 
The initial phase of the investigation involved the derivation of the nonequilibrium defect 
stream function formulation for incompressible flows. The formulation reduces to that of 
Mellor and Gibson for equilibrium flow. The primary accomplishment of the 
incompressible studies was, however, the development and implementation of the 
technique used to match the empirical inner region solutions with numerically computed 
outer region solutions. Another new feature was the use of a law of the wake in addition to 
the law of the wall. An existing law of the wake was used in this part of the investigation. 

The primary phase of the investigation was the development of the nonequilibrium, 
compressible defect stream function formulation and the corresponding compressible inner 
region treatment in conjunction with the matching of the inner and outer region solutions. 
The compressible formulation was designed to reduce to the incompressible form of the 
initial work. A method for determining the coefficients of a postulated compressible law of 
the wake was also developed. This resulted in an analytic equation for the compressible 
law of the wake valid from the wall through the inner part of the outer region. 

The final phase consisted of the application of the present techniques in primitive 
variable form. The method was incorporated into an existing two-dimensional Navier- 
Stokes code and tested for several cases of compressible flow over a flat plate. 

4.1 Incompressible Flow 

Solutions for incompressible, equilibrium boundary layers have been computed 
with the asymptotic, or zero-order, and full-equation forms of the present method. As 
discussed earlier, the zero-order form is taken in the limit of vanishing shear- stress velocity 
ratio. All solutions of this section were computed using the incompressible law of the 
wake of Moses which is one of the widely accepted incompressible laws along with that 
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of Coles 20 . The Moses law was chosen because of its simple, polynomial form. Moses' 
incompressible law of the wake is 

h(tn)=-n<p) [3(^-) 2 - 2(^) 3 ] 
k 5 6 

where 5 is the boundary layer thickness. The ratio (A/5) and the coefficient n are functions 

22 

of (3, and P is a function of the streamwise coordinate White gives empirical 
expressions for these terms with regard to incompressible, equilibrium turbulent flow as 
follows: 

A (1 +fl) 

6 K 


and 

n=^(o.5 + p) 3/4 (4.1) 

The law of the wall given in equation (2.8) reduces exactly to the common logarithmic-law 
form. The empirical coefficients for the computations of this section were chosen to 
conform with the work of Mellor and Gibson ; the law of the wall constants k and B have 
the values 0.41 and 4.9, respectively, and the arbitrary function K(x,y) in the outer region 
eddy-viscosity model is assumed to be the Clauser constant and has the value 0.016 after 
Mellor and Gibson. 

The equilibrium boundary layer problem is solved easily using a shooting 

technique. A value for the match point velocity defect fm is iterated until the far-field 

boundary condition on f (equation(2.16)) is satisfied. With the guess of fm, the match 

point location and then the corresponding value of f m are determined. A fourth order 

22 

Runge-Kutta integration routine given by White is used to integrate from the match point 
across the outer region of the boundary layer. 
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The present solutions are compared with the full-equation solutions of Mellor and 

32 

Gibson . All of these results are for a Reynolds number based on the edge velocity and 
displacement thickness (Reg*) of 10 5 . The comparison shown in figure 3 includes the 
zero-order closed form solution for the most favorable pressure gradient, P=-0.5. This 
solution is shown in Appendix B. The present zero-order numerical solution is seen to be 
in complete agreement with the closed form solution and in close agreement with Mellor 
and Gibson's full-equation solution. The zero-order solutions provide an indication of the 
accuracy of the numerical scheme; the numerical solutions give an early indication of the 
accuracy of the zero-order equations. 

It should be noted in figure 3 that the present method has resolved a turbulent 
boundary layer with only eleven grid points across the boundary layer. It is fair to question 
the accuracy of a solution with so few grid points. Figure 4 compares various grids for a 
typical value of P (P = 2) and shows that the eleven point grid solution is in agreement with 
the 51 and 101 point solutions. Only the six point solution deviates noticeably from the 
fine grid solutions, and this deviation is negligible. On the basis of these observations, the 
present solutions were computed with eleven uniformly spaced grid points. The first grid 
point is always at the match (patch) point, and the edge value of the normal coordinate is 
fixed. This is the adaptive grid approach mentioned in section 3.3. The normal coordinate 
reduces to Tj for incompressible flow and is plotted as such. During the iteration of the 
velocity defect f^, the first grid point fj m shifts with fm according to equation (3.6). The 
uniform grid is recomputed for each value of f^. This grid procedure was used for 
convenience rather than necessity as the grid points other than the first could have been 
fixed. By eliminating the inner region computations, a grid-point savings on the order of 
fifty percent was realized as was the case in references 12, 13, 17, and 18. 

Equilibrium velocity defect profiles are presented in figure 5 for several small-p 
cases (P < 1); and in figure 6, results are presented for large-P cases (P > 1). The 
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distinction between small and large (5 cases was made by Mellor and Gibson because the 
standard velocity defect variables become unbounded as P approaches infinity and the 
shear-stress velocity u* approaches zero. They also observed the edge value of the normal 
coordinate rj e also approaches zero as p approaches infinity. Mellor and Gibson presented 
a transformation to overcome these difficulties. The transformation effectively replaces the 
shear-stress velocity u* with a "pressure velocity" u p which is defined as 

u p = 

for incompressible flow. This quantity is related to u* and p as 
up = u* VF 

The incompressible coordinates are defined as 
S = Ps and N = rj Vp" 

and the velocity defect normalized by the pressure velocity and the transformed shear-stress 
velocity ratio are 

and X = WP 

3N U P t/|3 

respectively. Although it was not necessary to invoke this transformation in order to 
achieve the solutions presented in this section, the variables F 1 and N are used for the 
large-P cases to facilitate comparison with the results of Mellor and Gibson. 

It can be seen in these figures that the full-equation form of the present method 
compares extremely well with the full-equation solutions of Mellor and Gibson. Even 
more important is the excellent agreement of the present zero-order solutions with the full- 
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equation solutions. The relative simplicity of the zero-order approach and its excellent 
results makes it a desirable alternative to the full-equation approach. 

The differences seen in the inner region are caused by the differences between 
Mellor and Gibson's inner region eddy-viscosity model and the present combination of the 
law of the wall and law of the wake. Mellor and Gibson's solutions approach the law of 
the wall in the limit as fj goes to zero. As will be seen later, Mellor and Gibson's approach 
allows a large deviation from the law of the wall and wake in the inner region. This effect 
becomes more noticeable with increasingly strong adverse pressure gradients (increasing 

P). 

The P = 10 case shown in figure 6 was computed with the Clauser constant 
k = 0.0154935 rather than k = 0.016. As demonstrated in figure 7 for the zero-order 
approach, a match point does not exist for p = 10 with k = 0.016 because the far-field 
boundary condition Foo = - 1 cannot be satisfied. This failure to obtain a solution indicates 
the problem is overspecified. Figure 7 shows that the value k = 0.0154935 yields a single 
value of Fm that satisfies the far-field boundary condition for P = 10. The maximum value 
of P for which a solution exists for k = 0.016 is approximately 8.6. For values of P less 
than 8.6, there are two values of fm which allow the far-field boundary condition to be 
satisfied. The physically correct value is the smaller of the two. 

The problem of overspecification is alleviated by allowing the law of the wake to be 
computed as part of the solution, as discussed in section 4.2, rather than relying on a fixed 
law of the wake as done when computing the solutions shown in this section. The 
flexibility of computing the law of the wake allows solutions to be computed without 
altering the Clauser constant as was done for the P = 10 case shown in figure 4. 

Results for the shear-stress velocity ratio y from the present solutions are compared 
with the results of Mellor and Gibson in figure 8. The two sets of results are in close 
agreement but are not identical. 
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The present method for evaluating y differs substantially from that of Mellor and 
Gibson , Walker et al , and others in that the value of y is obtained from a patching of 
the inner and outer solutions at a particular point in the present method whereas a matching 
of the outer limit of the inner solution with the inner limit of the outer solution is used in the 
other methods. These methods are not strictly asymptotic in that the inner and outer eddy- 
viscosity models are patched at a point analogous to the patch point used in the present 
method. The equation Mellor and Gibson used to evaluate y is very similar to equation 
(2.1 1), the equation used in the present method. The difference is one of application; 
equation (2.1 1) is evaluated at a point a finite distance from the surface in the present 
method, while the equation used by Mellor and Gibson is evaluated as close to the surface 
as possible (the inner limit of the outer solution). As a result, the values of y obtained by 
Mellor and Gibson are consistently smaller than those obtained by the present method. 
However, larger values of y can be obtained from the results of Mellor and Gibson by 
evaluating their equation at the point where the eddy viscosities are patched. 

The inner region velocity profiles of the present method differ from those of Mellor 
and Gibson. The differences increase with an increasingly strong adverse pressure 
gradient (increasing (3) and are best displayed using inner variable coordinates. Figure 9a 
shows the {3 — 4 case. It is clear in this figure that, although the inner solutions are 
identical near the wall, Mellor and Gibson's solution deviates from the logarithmic behavior 
at least an order of magnitude (in y + units) closer to the wall than in the present method. 
This shows the difference between applying an eddy-viscosity model that approaches the 
logarithmic form when approaching the wall and enforcing the logarithmic behavior 
throughout the inner region. In view of this difference, a comparison with experimental 
data seems appropriate. Figure 9b compares the present solution with experimental data 
taken from reference 46 for Causer's second equilibrium flow. The data shown in 
figure 9b correspond to profile 2305 in reference 46. For this profile, f3 = 7.531 and 
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Reg* = 30692.5. Figure 9b shows excellent agreement between the experimental data the 
present solution. Given the relatively large value of (3, one would expect less agreement, 
especially in the logarithmic region, between a Mellor and Gibson solution and experiment. 
The present solution and experiment were also in close agreement for y where a 0.87 
percent difference was found. Present solutions were computed and compared with 
experimental data from reference 46 for many values of (3 and good agreement was 
generally observed. Table 1 shows several comparisons of the experimentally and 
numerically determined shear-stress velocity ratios where the numerical solutions were 
computed with the zero-order formulation. The analytic solution from the method of 
reference 42 (see Appendix A) is also shown. 

Another feature of the present method can be observed quite well in figure 9b. 
Notice that the match point of the present solution is well beyond the logarithmic region. 
This result demonstrates the advantage gained when including the law of the wake in the 
present formulation. For all present solutions, the law of the wake allows the match point 
to move beyond the logarithmic region thus extending the usefulness of the analytic inner 
region treatment and further reducing the region which must be resolved numerically. As 
one would expect, the influence of the law of the wake increased as (3 increased. 

4.2 Compressible Flow 

Solutions for compressible, equilibrium boundary layers have been computed with 
the constrained zero-order and full-equation forms of the present method. The full- 
equation solutions, however, were generated with the equations of the matching procedure 
truncated to second-order. As discussed earlier, the constrained zero-order form is not the 
strict asymptotic form taken in the limit of vanishing y as the term e, which is of order yM e , 
is not always a small parameter and thus remains in the formulation. The significance of 
the term e increases with increasing M e and has been included in all of the present 
computations for compressible flow. The term, however, is practically negligible until 
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reaching the high supersonic and hypersonic flow regimes. To fulfill the goal of 
determining a compressible law of the wake, solutions were computed with an iterative 
curve fit technique supplying the coefficient w of the postulated law of the wake given in 
equation (2.9). An analytic expression for the coefficient w is determined from these 
results. 

The compressible equilibrium boundary-layer problem is solved using a shooting 
technique in the same fashion as described for incompressible flow. The primary 
difference is the determination of the law of the wake as part of the computational 
procedure. The specific procedure for determining the coefficient w is discussed later, 
suffice it to say that a new value of the coefficient w is determined with each iteration of fm. 
The empirical parameters k and B for the law of the wall and the outer region eddy- 
viscosity parameter k retain their values from the incompressible flow cases; they are 0.41, 
4.9, and 0.016, respectively. Note that the additional empiricism of the Crocco 
relationship is also incorporated in the compressible flow computations. 

Both the constrained zero-order and full-equation solutions have been computed 

A 

over M e , p, and Res* ranges. The solutions shown were computed with 32 uniformly 
spaced grid points across the outer region (the first point being the match point) in order to 
facilitate the law of the wake curve fit procedure. Fewer points could have been used with 
the same success if the points had been nonuniformly spaced and clustered near the match 
point. Significantly fewer points were needed when the analytic law of the wake was used 
(ie once the compressible law of the wake was determined and the curve fit procedure no 
longer needed). The inner region solutions shown in the following figures were 
determined analytically from the laws of the wall and wake and the Crocco relationship. 

Figure 10 compares velocity defect profiles for the constrained zero-order and full- 
equation formulations over a Mach number range for Reg* = 10 4 and P = 0 (the flat-plate 
case). It can be seen that the agreement between the formulations is extremely good. Table 
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2a shows the shear-stress velocity ratio y computations to be in agreement on the order of 
one percent across the Me range. The percent error of the constrained zero-order 
formulation for the prediction of y, as compared to the full-equation formulation prediction, 
increases as the edge Mach number decreases, but we have previously seen that the 
predictions show excellent agreement for incompressible flow. 

Figure 1 1 shows the comparison over a Reynolds number (Reg*) range for M e = 3 

A 

and P = 0. Again, the velocity defect profiles are nearly identical and, as shown in Table 
2b, the shear-stress velocity ratios as computed by the constrained zero-order formulation 
show an error on the order of one percent. As the Reynolds number decreases, the percent 
error of the constrained zero-order prediction of y increases. 

A 

Figure 12 shows the comparison for several P values for M e = 3 and Reg* = 10 4 . 

A 

The P = 1 case represents a relatively strong adverse pressure gradient case for supersonic 

A 

flow. The velocity defect agreement for the P =1 case, although not as good as for the flat 
plate case, is still quite good. Table 2c shows the comparison of the shear-stress velocity 
ratios computed by the two formulations. The percent error of the constrained zero-order 
formulation increases with increasingly strong adverse pressure gradient. Although the 
error seen here is larger than for variations of Me and Reg*, the error over this practical 

A 

range of P's remains below five percent. 

The reduced complexity of the constrained zero-order formulation and the excellent 
agreement of the results above justify it as an excellent approximation and alternative to the 
full-equation formulation. All subsequent solutions shown in this section are based on the 
constrained zero-order formulation. 

The determination of the law of the wake coefficient w is an important part of the 
present solution procedure. The assumed, or postulated, functional form of the law of the 
wake is shown in equation (2.9), where w is the unknown coefficient that is determined by 
a least squares curve fit during each iteration of fm. Recall that for the present treatment the 
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law of the wake need only be valid in the inner most part of the outer region; the curve fit to 
determine w focuses on this region of the boundary layer. In effect, a value for the law of 
the wake h is determined at each point in the outer region as 


h(n) = 



pw df 1 

— rr-g(n)+- 

Pe dn y 


(4.2) 


where g is the law of the wall, y is the current value of the shear-stress velocity ratio, and 
f is determined from the numerical solution. A least squares curve fit of the values of h is 
made with the value of the coefficient w as a result. A standard least squares correlation 
factor was used to judge the accuracy of the fit. The number of points used for the curve fit 
varied from case to case in order to maintain approximately the same correlation factor, that 
is, approximately the same level of accuracy with regard to the fit was maintained from case 
to case. Note that this procedure can be used for any general postulated law of the wake 
with any number of unknown coefficients. More complex laws of the wake which could 
be valid over more of the outer region of the boundary layer can be determined. 

Figure 13 compares the resulting analytic law of the wake as defined by equation 
(2.9) and the determined value of the coefficient w with the numerically computed values of 
the law of the wake as determined by equation (4.2) for several Mach numbers. The 
profiles shown are for Reg* = 10 4 and P = 0. It can be seen that the number of points 
accurately fit by the postulated law of the wake decreases as the Mach number increases; 
the fit is always accurate near the match point. From this figure, as well as other solutions 
not presented, it appears that the y 2 dominance of the law of the wake in the inner region 
begins to diminish as Me approaches 4. 

It can be seen in figure 1 3 that the law of the wake coefficient w is a function of Me 

A 

in addition to being dependent on p. Figure 14a shows the variations of w as a function of 

A A 

Me and P for Reg* = 10 4 . It is seen that w increases significantly as P increases, especially 
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at small values of Mg. As the Mach number approaches 4, w becomes a relatively small 

A 

number for all P, and the strong dependence of h on r\ 2 in the inner part of the boundary 
diminishes. Figure 14b shows the variations of the coefficient w as a function of Me and 

A 

Reg* for P = 0; the coefficient w is a weak function of Reg*. 

An analytic model for the law of the wake coefficient based on the above solutions 
can be written as 

w(M e ,P) = w(0,P) (£w) a( ^ 

Pe 

A 

where w(0,P) is the law of the wake coefficient for incompressible flow and the exponent 

A 

a is a function of P only. An equation for the coefficient w in incompressible flow follows 
from section 4.1 as 

w(o.P) = n (|) 2 = n (4.3) 

5 K 

22 

where the coefficient n is a function of p. White's equation for n, as given in equation 
(4.1), can be used. Values of equation (4.3) based on equation (4.1) are shown in figure 
14a in comparison to the values determined by the curve fit procedure. The present zero- 
order solutions have been used to determine a new expression for n. The equation is 

n = <^(0.5 + P)° 2 (4.4) 

where and cl, are functions of P, rather than constants as in equation (4.1), and are 
defined as 

and ° 2 4 + §«p[> 5+ e> M ] 
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Figure 15a compares equations (4.1) and (4.4) with values of IT deduced from the 
numerically determined values of the coefficient w. Figure 15b compares the 
corresponding values of the coefficient w. Equation (4.4) is seen to be an excellent 
approximation to the numerical data throughout the range of p. White never intended 
equation (4.1) to be used for high P cases, so the improvement in this range is not 
unexpected. The agreement between equations (4.1) and (4.4) is much better in the low P 
range, but the new expression does fit the numerical data better. The equation for a was 
determined from the solutions shown in figure 14a and is written as 

a(P) = 2-exp[-|V2p+ 1 ] ( 4 - 5 ) 

The resulting curve is shown in figure 16. 

As mentioned previously, significant grid reduction can be realized when using the 
analytic law of the wake. Figure 17 shows this grid reduction when the analytic equation is 
used rather than the curve fit procedure. The solutions shown for the case where Me = 3, 

P = 0, and Reg* = 10 4 are virtually identical with only the seven point grid solution 
deviating ever so slightly. The computed values of y were also essentially identical as the 
difference between the 32-point and the 7-point solutions was less than a quarter of one 
percent. 

A 

Figure 18 presents a series of constrained zero-order solutions over a P range from 
-0.3 to 2 for Reg* = 10 4 and M e = 3. Profiles of the velocity defect, velocity ratio, and the 
density ratio are given. Since it is unnecessary to evaluate the inner region solution in order 
to compute the outer region solution, only the outer region results are shown. Figure 19 
presents a series of flat plate solutions for a Mach number range from 0 to 4 and for 
Res* = 10 4 . Figures 18 and 19 are presented to document the effect of parameter variations 
on the profiles. It is interesting to observe that Me variations have a minimal effect on the 
velocity defect profiles. The velocity defect profiles for variations of Reg* with fixed Me 
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and P, although not presented here, were virtually identical. The fact that the velocity 
defect profiles are relatively insensitive to variations of Me and Res* supports the definition 

A A 

of P as the compressible equilibrium flow parameter. The parameter P effectively 
designates the various families of similar solutions for compressible equilibrium flow, at 
least under adiabatic wall conditions, in the manner shown by Clauser with his parameter P 

A 

for incompressible flow. Note that the present parameter p reduces to Clauser's P for 
incompressible flow. 

A 

Figure 20a shows the effect of Me and P on the shear-stress velocity ratio for 
Reg* = 10 4 . It is observed that y increases with increasing Me while decreasing with 
increasing P (increasingly strong adverse pressure gradient). For any combination of Me 
and p, y increases as the Reynolds number decreases. Figure 20b demonstrates this trend 
for flat plate flow. The trend was documented in reference 32 for incompressible flow. 

A 

Figure 21 shows the effect of Me, P, and Reg* on the compressibility parameter co. 
Recall that this parameter is the ratio of the velocity thickness to the density-weighted 
velocity thickness and appears throughout the present formulation. Figure 21a shows the 
Mach number and pressure gradient effects on to for Reg* = 10 4 . By definition, CO remains 
one for all incompressible flows. It is observed that <0 increases with increasing Me as well 

A 

as with increasing p. The effect of the pressure gradient is magnified with increasing Me. 

A 

Figure 21b shows the effect of Reg* on co for Me = 3 and P = 0. It is observed that co 
decreases with increasing Reynolds number. The curves of figure 21a would shift 
upwards for Reg* less than 10 4 and downwards for Reg* greater than 10 4 with the 
exception of the incompressible flow curve which remains fixed. 

Comparisons of the constrained zero-order solutions with experimental data 
compiled by Femholz and Finley 4 ^ have also been made. Figure 22 shows the profile 
comparisons for two flat plate cases; each numerical solution was computed with 32 grid 
points across the outer region and the curve fit procedure was used to determine the law of 
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the wake coefficient. The agreement of the velocity defect and velocity ratio profiles is 
good, while the density ratio profiles are in reasonable agreement. The use of the Crocco 
relationship for the density ratio magnifies the velocity differences, especially near the edge 
of the boundary layer. The modification of the outer eddy-viscosity model by including the 
intermittency factor would improve the comparisons because the gradient of f' would 
increase as the boundary layer edge is approached. An intermittency factor could be 
incorporated into the formulation easily by allowing the arbitrary function K in equation 
(2.4) to be variable rather than using the Clauser constant k. 

Table 3 compares the computed and experimental shear-stress velocity ratios. The 
computed solutions are from the constrained zero-order formulation and the experimental 
values are based on the values of skin friction as measured by a force balance. The analytic 
solutions shown were determined with the method of reference 42. The agreement for the 
flat plate cases is excellent; the experimental and computed values differ by no more than 
two percent. The 5801 series of experimental data represents a nonequilibrium flow in 
which data was taken in both favorable and adverse pressure gradients. The agreement 
shown for this series is not as good as seen with the flat plate cases. However, the largest 
difference occurs with the larger of the two adverse pressure gradient cases and the 
difference approximately 6.5 percent. The editors of reference 47 cast some doubt on the 
accuracy of the experimental values due to a lack of upstream history information and 
unknown conditions at the edge of the boundary layer. The value of the shear-stress 
velocity ratio should decrease between profiles 5801/0501 and 5801/0801 as Me decreases, 

A 

P increases, and Reg* increases; these trends, as observed in the present investigation, 
indicate a decrease in y that apparently was not observed experimentally. 

4.3 Primitive Variable Application 

The ultimate objective of the present investigation is implementation of the present 
concepts into existing three-dimensional, Navier-Stokes methods. It should be at this level 
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that the greatest rewards with regard to computational efficiency and accuracy will be 
realized. Coarse grids are often used near the wall in the interest of computational 
efficiency; the result is a loss of accuracy. Accurate numerical computations near the wall 
require a fine grid, which results in a loss of computational efficiency. The numerical 
computations near the wall can be replaced by an empirical treatment that does not rely on a 
fine grid for accuracy. The initial effort towards the final goal is discussed in this section. 

The concepts developed in the present investigation have been implemented into the 
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computer code known as NASCRIN (Numerical Analysis of Scramjet Inlets) ’ . The 
code solves the two-dimensional Navier-Stokes equations in conservation law form. 

Turbulence is modeled with the zero-equation eddy-viscosity model of Baldwin and 

in q 

Lomax . The zero-equation model of Cebeci and Smith has also been incorporated into 

NASCRIN and is used in the present investigation. A numerical coordinate transformation 

is used which generates a set of boundary-fitted curvilinear coordinates in order to facilitate 

the resolution of flowfields with general geometries. The transformation allows for the 

clustering of grid points in regions of high gradients. The transformed governing 

equations are solved using an unsplit, explicit, predictor-corrector, finite-difference method 

developed by MacCormack 5 ^. This finite-difference method is second-order accurate in 

both time and space. The flowfield variables at the inflow boundary are fixed at some 

designated initial condition; the variables at the outflow boundary are determined by a first- 

order extrapolation. No-slip and adiabatic wall conditions are applied on solid surfaces. 

The adiabatic wall condition used in NASCRIN @T/dy = 0 at y =0) has been replaced by a 

more accurate 51 condition (3H/5y = 0 at y =0) in the present investigation. 

The present inner region and matching treatments were incorporated into 
NASCRIN with relatively minor modifications to the code. The no-slip condition at the 
wall was replaced with the condition on the gradient of the velocity du/dy given in equation 
(3.7), and an additional subroutine was incorporated to determine the inner/outer match 
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point location, and the shear-stress velocity. A general description of the use of this 
boundary condition was given in section 3.3. The finite-difference method is used 
throughout the flowfield; only the outer region eddy viscosity of the zero-equation 
turbulence model (the Cebeci-Smith model, for example) is required. In essence, the 
resulting numerical solution is not physically accurate near the wall, as it is based on the 
outer region eddy-viscosity model and the velocity gradient boundary condition at the wall. 
However, the solution away from the wall is accurate, where the outer region eddy- 
viscosity model is physically correct. The physically accurate solution near the wall is 
determined from the empirical expressions used in the inner region treatment after the 
numerical solution of the outer region has converged. 

The velocity gradient boundary condition (equation (3.7)) depends on the value the 
shear-stress velocity u*. This parameter varies in the streamwise direction and is updated 
during the computation after every twenty time-steps. This updating scheme was chosen 
because NASCRIN updates the eddy-viscosity after every twenty time-steps. The shear- 
stress velocity is determined from the matching of the inner and outer region treatments. In 
this initial effort, the present concepts have been incorporated in the form of the method of 

A 

reference 42 (see Appendix A). The numerical method provides the values of Mg, (3, and 
Reg* needed in the matching subroutine, which is the new subroutine that has been 
incorporated into NASCRIN. This subroutine is, in fact, the method of reference 42, 
which solves analytically for the value of u* at each streamwise location. Note that the 
analytic solution of reference 42 is for the constrained zero-order equation with equilibrium 
flow. 

The NASCRIN code was used with and without the present modifications. 
Solutions with the present modifications are referred to as the slip-velocity solutions. 
Solutions computed without modification are referred to as the baseline solutions. The 
Cebeci-Smith turbulence model was used with the Clauser constant k = 0.0168; the law of 
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the wall constants K and B have the values 0.41 and 5.0, respectively, and are used in the 
present inner region and matching treatments. Baseline and slip-velocity solutions are 
compared with experimental data. 

The test cases considered in this initial effort were for compressible flow on flat 
plates. The two cases presented here are for M e = 2.578 and 4.544; test conditions and 
experimental results for these two flows are found in reference 47 as part of the 5301 series 
of data. These two cases were discussed previously in section 4.2. Inflow conditions 
were determined by the analytic method of reference 42 for a specified Reg*; the 
Me = 2.578 case was started with Reg* = 2500, and the M e = 4.544 case was started with 
Reg* = 400 at the inflow boundary. In each case, the experimental measurement station, 
corresponding to Reg* = 5295 and 902 for Me = 2.578 and 4.544, respectively, is near the 
middle of the computational field. The numerical solution was considered to be converged 
when one of the following criteria was met: (1) the maximum change in the density at each 
grid point was less than or equal to 10 -7 , or (2) the average change in the density over the 
entire field was less than or equal to 10' 8 . Each case was computed with fifty-one 
uniformly-spaced grid points in the streamwise direction and a constant Ax of 0.0168 
meters. 

The grid reduction, and corresponding increase in computational efficiency, is 
directly related to the grid-spacing normal to the surface. It was found that fifty-one grid 
points normal to the surface was adequate to resolve the flowfields for each case with the 
baseline method. The first grid point off of the wall was located at y = 0.000027 meters, 
which corresponds to y + values of approximately 3.9 and 0.9 for the M e = 2.578 and 
4.544 cases, respectively. Slip-velocity solutions were determined with far fewer grid 
points normal to the surface. To demonstrate, each case was computed with twenty-six 
points. The first grid point off of the wall is located at y = .000142 meters, or 
approximately five times further off of the wall than was used in the baseline solutions. 
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The values of y + were approximately 21 and 4.9 for the Me = 2.578 and 4.544 cases, 
respectively. The match point was located at y = 0.0021 and 0.00145 meters, with 
corresponding y^ values of 312 and 50, for the Me = 2.578 and 4.544 cases, respectively. 
Recall that from the match point outward, the slip-velocity solution is physically accurate; 
empirical expressions provide the physically accurate solution between the wall and the 
match point. 

Velocity ratio and density ratio profiles for the two cases are shown in figures 23 
and 24. The baseline and slip-velocity solutions show excellent agreement in the outer 
region; the empirical inner region solution, based on the matching process, is in excellent 
agreement with the inner region solution of the baseline method. The only significant 
deviation from experiment occurs near the edge of the boundary layer. The predicted 
values of the shear-stress velocity ratio for Mg = 2.578 were 0.0430 and 0.0423 for the 
baseline and slip-velocity solutions, respectively; the shear-stress velocity ratios for 
Me = 4.544 were 0.0522 and 0.0532 for the baseline and slip-velocity solutions, 
respectively. The corresponding experimental values were 0.0426 and 0.0544 for Me = 
2.578 and 4.544, respectively. In each case, the slip-velocity solution provides the better 
prediction. 

The time-step is directly proportional to the minimum grid spacing normal to the 
surface. As discussed earlier, the slip-velocity procedure enables an increase in the 
minimum spacing (by moving the first grid point off of the wall further from the wall). 

The effect on the rate of convergence is dramatic. The baseline solution for the Mg = 2.578 
case required 26250 steps, at an average At of approximately 0.00000008 seconds; the 
corresponding slip-velocity solution required only 5900 steps at an average At of 
approximately 0.0000004 seconds. The average time- step increases by approximately five 
for the slip-velocity solution. The baseline solution for the Me = 4.544 case required 
21250 steps, at an average At of approximately 0.00000008 seconds; the corresponding 
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slip-velocity solution required only 9000 steps at an average At of approximately 
0.0000004 seconds. Again, the average time-step increases by approximately five for the 
slip-velocity solution. 

Slip-velocity solutions were computed for a variety of grids; the twenty-six point 
solutions just discussed should not be taken as the optimum result. In fact, the slip- 
velocity solutions were relatively insensitive to the number of grid points normal to the 
surface. This statement holds for the resulting parameter profiles as well as the predicted 
shear-stress velocity ratio. The effect of reducing the amount of grid clustering (near the 
wall) and the number of grid points is to increase the rate of convergence. 
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5 Concluding Remarks 


Equations governing two-dimensional, nonequilibrium, compressible turbulent 
boundary-layer flow have been derived using a defect stream function formulation. This 
formulation is of particular interest because it has a constrained zero-order approximation 
for which the tangential momentum equation has a first integral. The approximation is 
obtained in the limit of vanishing shear-stress velocity to edge velocity ratio and the 
constraint exists because the parameter yMe cannot be neglected for high-speed flows. The 
modified Crocco temperature-velocity relationship has been used as an approximation to the 
energy equation. Compressible equilibrium flow has been defined with this formulation. 
Solutions for incompressible and compressible equilibrium flows have shown the 
constrained zero-order form to compare well with the full-equation form and experimental 
data. 

The lack of popularity of previous defect stream function formulations is probably 
due to the difficulty encountered in enforcing the no-slip surface boundary condition. The 
present method overcomes this difficulty with a law of the wall/law of the wake 
formulation for the inner part of the boundary layer. This inner region treatment is 
mathematically patched to the outer region formulation. This formulation eliminates the 
need for an inner region eddy- viscosity model. The method also significantly reduces the 
grid requirements by eliminating numerical computations in the inner region. 

A compressible law of the wake has been developed. The law of the wake was 
assumed to be proportional to y 2 and an iterative curve fitting procedure was used to 
determine its coefficient as part of the computation. Solutions generated with this technique 
allowed the definition of an analytic expression for the coefficient, and thus a completely 
analytic compressible law of the wake. Results show that use of the law of the wake 
permits patch point locations beyond the logarithmic region of the boundary layer. 


58 



The method, in its constrained zero-order equilibrium flow form, has been 
implemented into an existing two-dimensional Navier-Stokes code. Solutions for 
compressible, turbulent boundary-layer flow on flat plates have been computed. As 
implemented, solutions were computed using only an outer region eddy- viscosity model 
from the surface to the edge of the boundary layer. The no-slip boundary condition for the 
velocity at the surface was replaced by a boundary condition on the normal derivative of the 
tangential velocity at the surface. Solutions showed a significant increase in computational 
efficiency as a result of a coarse grid near the wall and increased time steps. The solutions 
showed excellent agreement with experimental data, analytic solutions, and the baseline, 
Navier-Stokes computations. 

Future investigations should include the incorporation of the nonequilibrium, defect 
stream function equations and matching treatments into NASCRIN, in addition to other 
Navier-Stokes solvers with various turbulence models. The current defect stream function 
treatment should be modified to account for nonadiabatic wall effects. Finally, the 
treatment should be extended to and implemented in three dimensions. An initial look at the 
extension to three dimensions, for adiabatic wall conditions, is given in reference 52. 
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6 Appendices 


6.1 Appendix A: Analytic Solution for Compressible Turbulent 

Flow 

The constrained zero-order analytic solution for compressible, turbulent boundary 
layers and resulting skin friction law are discussed in detail in reference 42. This appendix 
provides a general outline of the solution. 

As discussed in section 2.7, the constrained zero-order equation for equilibrium 
flow given in equation (2.21) is essentially linear from incompressible to supersonic flow, 
but an additional transformation is necessary to obtain a linear governing equation for 
hypersonic flow. The transformation is 

s = s and fj= f — dr| = rj - ef 

J p 

The transformed governing equation is 

^=coK^+(l+2P)(fl-f| 0 ) ^.-f-i- e pf|L (6.1) 

ds aXy- df| drj 

where f is understood to be f 0 and 

P = 2p-ii^Ji(££--l) and = 

Pw (l+2p) 

This equation is linear if the inequality 

e Ipl « 1 (6.2) 

is satisfied. This inequality is satisfied for all pressure gradients from incompressible to 
supersonic flow and for small pressure gradients in hypersonic flow. An equilibrium form 
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of equation (6.1) can be determined in the same manner as discussed in section 2.8 with the 
additional condition that the inequality of equation (6.2) is satisfied. The equilibrium form 
of equation (6.1) is 

cok^ + (1+2P) (fl - %) — - f - 1 = 0 (6.3) 

dfj 7 dr) 

where the arbitrary function K is assumed to be the Clauser constant k. 

32 

As mentioned in section 2.8, Mellor and Gibson showed that the incompressible 
equation is in the form of a confluent hypergeometric equation. Following Mellor and 
Gibson, the present independent variable is defined as 

XT 1+2$ (tS - tTq) 2 

N s 

G)k 2 


The dependent variable for the present problem is defined as 
f=(l +f)exp(N) 

which has the free-stream boundary condition 

A 

fco=0 


In terms of these variables, equation (6.3) is 


XT d 2 f , ,1 XTX df 1 1 \* n 

N JX ,9 + (9 - N ) hn " 9 A 1 + A/f-0 


dN 2 


dN 2 


1 + 2p 


Confluent hypergeometric equations are of the form 


43 


-A A 

XT d 2 f . M . df 
N — r+ (b - N) ^rr 
dN 2 dN 


af = 0 


The parameters a and b for the present problem are 
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1 


b= I 


and 


4o 


-) 


1 + 2(5 


and the solution is 

f = C M(a,b,N) + D N 1 ' 5 M(l+a-b,2-b,N) 


where C and D are coefficients related to the boundary conditions and M(a,b,N) is 
Rummer's function which is given as 43 


M(a,b,N) = 1 + + 


a(a+l) N 2 
b(b+l) TT 


+ ... 


The resulting equation for the velocity profile in the outer region of the compressible, 
equilibrium turbulent boundary layer is 


^ - C { 3 M(a - 1/2, 1/2.N) - 

df l- V( a - 1/2) tok r ( 1+a ) 


— M(a,3/2,N) } exp(-N) 

cok 


where T is the gamma function and the coefficient C is given in reference 42. 

The inner region empirical formulation is identical to that previously described. The 
matching of the velocity and velocity gradients of the inner and outer region solutions gives 
an equation for the match point fj m , and finally an equation for the shear-stress velocity 
which is 


{£- { -l n (^\/^ fw) + B + - W — 

“e k mV n... K p w 


+ ^._Pe_ /£ w £ a r(l/2 + a) 

Pm \ Pe V(a- l/2)cok r (! + a ) 

(fjm ■ ffn) I i-I 

— ^M(a,3/2,N m ) ] exp(-N m )j 

cok 


M(a - l/2,l/2,N m ) 
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This equation is written in terms of the simplified form of equation (2.8) for the law of the 
wall. Solutions presented here, however, were computed using the full van Driest form of 
the law of the wall as given in equation (2.8). 


6.2 Appendix B: Analytic Solution for Incompressible 

Turbulent Flow at P =-0.5 

A closed-form solution exists for the zero-order equilibrium flow case with regard 
to the most favorable pressure gradient (P = -0.5). The zero-order governing equation for 
equilibrium flow with P = -0.5 and K defined as the Clauser constant k is 

kfo-f 0 =l 

Evaluation of equation (4. 1) for n shows that the law of the wake vanishes for P = -0.5. 
Therefore, equation (3.2) has the form 

f 0,m = ; nm 0o,m ) 

K. 

The equation for tj m is 

= Vl +‘*f 0m k/(fi m K) - 1 

Mm — 

2f 0,n/ *Tm 

The solution is 

f 0 Olm)=— exp[(ri m -fi)/Vk] - 1 

Kfim 

where 

Tim = 0.5 [ K - Vk - V(K - Vk ) 2 - 4k ] 
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The velocity defect is 


^=6 = — 

U Kn m 

and the shear-stress velocity is 


^ [-In (Reg* ilm) + B + 

K 


- fj) / Vk ] 


Vk j 

KTjm 
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Case 1 " 

Me 

A 

P 

Res* y 

experiment 

Y 

numerical 

Y 

Jfcsfc 

analytical 

1312 

0 

-0.230 

5447.0 

.0411 

.0407 

.0401 

1410 

0 

0 

6754.2 

.0388 

.0386 

.0382 

1422 

0 

0 

19,321.1 

.0350 

.0351 

.0348 

2205 

0 

1.891 

20,230.2 

.0319 

.0301 

.0296 

2305 

0 

7.531 

30,692.5 

.0231 

.0229 

.0225 

2402 

0 

3.026 

42,832.9 

.0267 

.0265 

.0265 


* experimental data from reference 46. 

analytical solution from method of reference 42. 

Table 1. Shear-stress velocity ratio comparisons with experimental data 
and analytical solutions for incompressible flow. 
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Me 

A 

P 

Reg* 

Yfull eqn 

^zero-order 

% error 

0 

0 

10 4 

.0367 

.0372 

1.362 

1 

0 

10 4 

.0375 

.0380 

1.333 

2 

0 

10 4 

.0392 

.0397 

1.276 

3 

0 

10 4 

.0409 

.0414 

1.222 

4 

0 

10 4 

.0424 

.0428 

0.943 

5 

0 

10 4 

.0435 

.0439 

0.920 


a) The effect of Me variations. 


Me 

A 

P 

Reg* 

Yfull eqn 

^zero-order 

% error 

3 

0 

10 3 

.0511 

.0523 

2.348 

3 

0 

10 4 

.0409 

.0414 

1.222 

3 

0 

105 

.0341 

.0344 

0.880 

3 

0 

10 6 

.0293 

.0295 

0.683 

3 

0 

10 7 

.0256 

.0258 

0.781 

b) The effect of Reg* variations. 

Me 

A 

P 

Reg* 

Y full eqn 

^zero-order 

% error 

3 

-0.15 

10 4 

.0418 

.0421 

0.718 

3 

0 

10 4 

.0409 

.0414 

1.222 

3 

0.5 

10 4 

.0386 

.0394 

2.073 

3 

1 

10 4 

.0367 

.0378 

2.997 

3 

2 

10 4 

.0339 

.0352 

3.835 


A 

c) The effect of (3 variations. 

Table 2. Shear-stress velocity ratio comparisons of constrained zero-order 
and full equation equilibrium solutions. 
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Case* 

Me 

A 

P 

Res* 

Y 

experiment 

Y 

numerical 

Y 

analytical 

5501/0101 

1.724 

0 

4719 

.0414 

.0420 

.0414 

5301/0601 

2.578 

0 

5295 

.0426 

.0431 

.0424 

5301/1302 

4.544 

0 

902 

.0544 

.0555 

.0533 

5801/0101 

1.947 

-0.155 

7949 

.0397 

.0412 

.0407 

5801/0501 

1.918 

.2175 

8596 

.0395 

.0391 

.0384 

5801/0801 

1.873 

.7265 

9024 

.0395 

.0370 

.0360 


* experimental data from reference 47. 

** analytical solution from the method of reference 42. 

Table 3. Shear-stress velocity ratio comparisons with experimental data 
and analytical solutions for compressible flow. 
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Figure 1. Laminar and turbulent velocity profiles on a flat plate. 

(reprint from reference 1 1) 
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a) Several experimental turbulent velocity profiles on flat plates. 
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b) Universal plot of experimental turbulent velocity profiles on flat plates. 

figure 2. Laminar-like behavior of the outer region of turbulent boundary layers. 

(composite of reprints from reference 3 1 ) 
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Figure iL.. Grid resolution comparison for incompressible flow. 
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a) (1 = 4/3 

Figure 6. Velocity defect profiles with (1 > 1 for incompressible flow. 
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b) p = 2 

Figure 6. Continued. 
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Figure 12. Effect of fi on the velocity defect profiles with M c = 3 and Reg* = 10 4 . 
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b) Effect of Reg* and M c for flat plate flow. 
Figure 14. Concluded. 
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a) Comparison of analytical and numerical values of n. 
Figure 15. The incompressible law of the wake coefficient. 
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based on eqn (4.4) 



b) Comparison of analytical and numerical values of w. 
Figure 15. Concluded. 
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b) Effect on velocity profiles. 
Figure 18. Continued. 
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c) Effect on density profiles. 
Figure 18. Concluded. 
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b) Effect on velocity profiles. 
Figure 19. Continued. 
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nsity profiles. 


Concluded. 
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a) Effect of p and M e with Res* = 10 4 . 

Figure 2 1 . Trends of the compressibility parameter co. 
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a) Comparison of velocity defect profiles. 

Figure 22. Comparison of various profiles with experimental data for compressible flow 

on a flat plate . 
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b) Comparison of velocity profiles. 
Figure 22. Continued. 
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c) Comparison of density profiles. 
Figure 22. Concluded. 
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a) Velocity ratio comparison. 

Figure 23. Baseline and slip-velocity Navier-Stokes solutions for M e = 2.578. 
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a) Velocity ratio comparison. 

Figure 24. Baseline and slip-velocity Navier-Stokes solutions for M e = 4.544. 
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b) Density ratio comparison. 
Figure 24. Concluded. 
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